%!TeX encoding = utf-8
%!TeX program = xelatex
\documentclass[12pt,a4paper]{article}                                                                            
\usepackage{amsmath,amsfonts,amsthm,amssymb}%AMS数学符号字体等的宏包
\usepackage{amsrefs}                           %ams参考文献宏包, 不需要的用户请删除此行!
%\usepackage{hyperref}                          %生成超链接
%\usepackage[colorlinks,linkcolor=black,anchorcolor=black,citecolor=black]{hyperref}

\usepackage{geometry}                        %页面设置
\usepackage{xeCJK}                               %中文字体宏包
\usepackage{graphicx}                          %插图
\usepackage{setspace}  
\usepackage{caption}
\captionsetup{font={small}}
\usepackage{listings} %C, C++等代码高亮，需要配合xcolor宏包支持
\usepackage{tabularx}
\usepackage{xcolor}

%%%%%%%%%%实用功能%%%%%%%%%%%%%%%%%%%%%%
\usepackage{indentfirst} % 解决 CJK 手段缩进问题,\setlength{\parindent}{2em}

%%%%%%%%%算法伪代码%%%%%%%%%%
\usepackage{algorithm}
\usepackage{algorithmicx}
\usepackage{algpseudocode}
\floatname{algorithm}{算法}
\renewcommand{\algorithmicrequire}{\textbf{输入:}}
\renewcommand{\algorithmicensure}{\textbf{输出:}}


\usepackage{cite}

%%%%%%%%%页面设置%%%%%%%%%%%
\setlength{\textwidth}{162mm}
\setlength{\textheight}{231mm}
\setlength{\headheight}{0cm}
\setlength{\topmargin}{-0.1cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{0cm}
\setlength{\parskip}{1mm}
\setlength{\unitlength}{1mm}
\setlength{\parindent}{2.06em}
\linespread{1.5}

%%%%%%%%%%重命名%%%%%%%%%%
\renewcommand\refname{参考文献}
\renewcommand\contentsname{目录}
\renewcommand{\figurename}{图}
\renewcommand{\tablename}{表}
\renewcommand\abstractname{\bf\large 摘~~要}


%%%%%%%%%%定义编号规则%%%%%%%%%%
\renewcommand{\thefigure}{\arabic{section}.\arabic{figure}}
\renewcommand{\thetable}{\arabic{section}.\arabic{table}}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}


%%%%%%%%%%CJK字体设置%%%%%%%%%%
\setCJKmainfont[BoldFont=STSongti-SC-Bold, ItalicFont=STKaiti]{STSong}
\setCJKsansfont[BoldFont=STHeiti]{STXihei}
\setCJKmonofont{STFangsong}

\newcommand\dd{\mathrm{d}}

%%%%%%%%%修改行间距%%%%%%%%%%%%
\usepackage{setspace}%使用间距宏包

\begin{document}

% Listing环境全局设置
\lstset{
%numbers=left,
%numberstyle=\tiny,
keywordstyle=\color{blue!70}, 
commentstyle=\color{red!50!green!50!blue!50},
frame=shadowbox,
rulesepcolor=\color{red!20!green!20!blue!20},
escapeinside=``,
xleftmargin=2em,
xrightmargin=2em,
aboveskip=1em
}





\def\chntoday{\the\year~年~\the\month~月~\the\day~日}

\title{并行求解布朗斯基方程}
\author{姚伟鹏\footnote{电话：13370130802，邮箱：ywp@pku.edu.cn，学号：1501110118 }}
\date{\chntoday}

\maketitle                            % 生成标题

\tableofcontents                      % 插入目录

\newpage

\section{问题模型}

\subsection{数学描述}

布朗斯基方程的原始变量形式为：		
\begin{equation*}
	\begin{split}
		&\frac{\partial \rho}{\partial t} + \vec{u} \cdot \nabla \rho = 0 \\
		&\frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} + \nabla \rho = {0 \choose \rho} \\
		&\nabla \cdot \vec{u} = 0
	\end{split}
\end{equation*}

其流函数—涡量形式为：
\begin{equation*}
	\begin{split}
		&\frac{\partial \rho}{\partial t} + \vec{u} \cdot \nabla \rho = 0 \\
		&\frac{\partial \omega}{\partial t} + \vec{u} \cdot \nabla \omega = -\frac{\partial \rho}{\partial x} \\
		&-\Delta \psi = \omega \\
		&\omega = \frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} \\
		&u = -\psi_y, v = \psi_x
	\end{split}
\end{equation*}

\subsection{求解要求}

\begin{itemize}
	\item 在正方形区域上用傅里叶拟谱方法求解布朗斯基方程；
	\item 在时间方向上使用高阶数值积分算子；
	\item 使用FFTW一维串行函数；
	\item 手工实现二维FFT，然后并行化
	\item 计算规模最低要求$512 \times 512$，期望做到$4096 \times 4096$。
	\item $t \in [0,3.16]$的初值为：
		\begin{equation*}
			\begin{split}
				&\omega(x,y;0) = 0 \\
				&\rho(x,y;0) = 50\rho_1(x,y) \rho_2(x,y) (1-\rho_1(x,y))
			\end{split}
		\end{equation*}
		其中：
		\begin{equation*}
			\rho_1(x,y) = \left\{
			\begin{array}{ll}
				\exp \left \{ 1-\frac{\pi^2}{\pi^2-x^2-(y-\pi)^2} \right \}. & \quad if \quad x^2+(y-\pi)^2 < \pi^2 \\
				0 & \quad otherwise.
			\end{array}
			\right.
		\end{equation*}
		\begin{equation*}
			\quad\rho_2(x,y) = \left\{
			\begin{array}{ll}
				\exp \left \{ 1-\frac{(1.95\pi)^2}{(1.95\pi)^2-(x-2\pi)^2} \right \}. & if \quad (x-2\pi)^2 < (1.95\pi)^2 \\
				0 & otherwise.
			\end{array}
			\right.
		\end{equation*}
\end{itemize}
	
%\subsection{物理预期}

\subsection{求解策略}
由于上述问题模型的原始变量形式和流函数—涡量形式可以证明是相互等价的，
而流函数—涡量形式更便于进行傅里叶拟谱方法。
所以我们考虑对流函数—涡量形式的布朗斯基方程进行数值求解。

在时间上，利用高精度的Runge-Kutta单步法对布朗斯基方程进行计算；
在空间上，利用傅里叶拟谱方法对布朗斯基方程进行计算。

\section{在空间方向利用傅里叶拟谱方法离散布朗斯基方程}

\subsection{傅里叶拟谱方法}

谱方法是一种重要的分析与计算工具。
在微分方程中，谱方法常可以将微分方程变换为代数方程，从而给计算带来方便。

谱方法的要点是把解展开成光滑函数（一般是正交多项式）的有限级数展开式。
根据基函数的不同，有不同的展开方式。
在这里，由于考虑的是周期性问题，所以在空间上选择以三角函数作为基函数的拟谱方法。

\subsection{流函数—涡量形式的Fourier变换}
将$\rho, \omega, u, v$按照如下方式进行Fourier展开：
\begin{equation}
	f(x_{j1},x_{j2};t) = \sum\limits_{k_1,k_2=-N/2}^{N/2-1}\hat{f}^t_{k_1 k_2}\exp{[i(k_1 x_{j1} + k_2 x_{j2})]}, \quad j_1,j_2 = 0,1,\cdots, N-1.
\end{equation}
其中，Fourier系数为：
\begin{equation}
	\hat{f}^t_{k_1,k_2} = \frac{1}{N^2} \sum\limits_{j_1,j_2=0}^{N-1} f(x_{j1},x_{j2};t)\exp{[-i(k_1 x_{j_1} + k_2 x_{j_2})]}, \quad k_1, k_2 = -\frac{N}{2}, \cdots, \frac{N}{2}-1. 
\end{equation}

记$\Phi = \vec{u} \cdot \nabla \rho$, $\Psi = \vec{u} \cdot \nabla \omega$，
由于Fourier变换具有如下的性质：
\begin{equation}
	\widehat{\nabla f} = i \vec{k} \hat{f}, \quad \hat{f} \otimes \hat{g} = \widehat{fg}
\end{equation}

所以我们可以将流函数—涡量形式写成如下方程组：
\begin{equation}
	\label{boussinesq_fft}
	\begin{split}
		& \frac{\partial \hat{\rho}^t_{k_1,k_2}}{\partial t} + \hat{\Phi}^t_{k_1,k_2} = 0 \\
		& \frac{\partial \hat{\omega}^t_{k_1,k_2}}{\partial t} + \hat{\Psi}^t_{k_1,k_2} = -i k_1 \hat{\rho}^t_{k_1,k_2} \\
		& (k_1^2 + k_2^2) \hat{\psi}^t_{k_1,k_2} = \hat{\omega}^t_{k_1,k_2} \\ 
		& \frac{\partial \hat{\omega}^t_{k_1,k_2}}{\partial t} = -i k_2 \hat{u}^t_{k_1,k_2} + i k_1 \hat{v}^t_{k_1,k_2} \\
		& \hat{u}^t_{k_1,k_2} = -i k_2 \hat{\psi}^t_{k_1,k_2}, \quad \hat{v}^t_{k_1,k_2} = i k_1 \hat{\psi}^t_{k_1,k_2}.
	\end{split}
\end{equation}

\subsection{利用FFT算法计算卷积}

在计算$\hat{\Phi}^t_{k_1.k_2}$和$\hat{\Psi}^t_{k_1,k_2}$时，需要计算卷积。
若直接在谱空间内进行计算，将会耗费非常大的计算量。
为了避免这一问题，我们通过离散Fourier逆变换，得到$\vec{u}$和$\nabla \rho$以及$\nabla \omega$在实空间的值。
这样一来，我们直接在实空间用乘法计算$\vec{u} \cdot \nabla \rho$以及$\vec{u} \cdot \nabla \omega$的值，
然后再通过离散Fourier变换将其结果投影到谱空间，从而得到$\hat{\Phi}^t_{k_1.k_2}$和$\hat{\Psi}^t_{k_1,k_2}$的值。

这种利用FFT算法来计算卷积的方法可以极大的减少计算量。

\subsection{手工实现二维FFT并行}

由于在第一章第二节中，要求使用FFTW一维串行函数来手工实现二维的FFT并行，
所以首先在一个方向上进行离散Fourier变换得到中间值$\hat{f}^t_{k_1,j_2}$或者$\hat{f}^t_{j_1,k_2}$。
然后再在另一个方向再进行一次离散Fourier变换得到$\hat{f}^t_{k_1,k_2}$。如下式所示：
\begin{equation}
	\begin{split}
		& \hat{f}^t_{k_1,j_2} = \frac{1}{N} \sum\limits_{j_1=0}^{N-1} f(x_{j_1},x_{j_2};t) \exp{(-i k_1 x_{j_1})}  \\
		& \hat{f}^t_{k_1,k_2}= \frac{1}{N} \sum\limits_{j_2=0}^{N-1} f(x_{k_1},x_{j_2};t) \exp{(-i k_2 x_{j_2})}
	\end{split}
\end{equation}

\begin{equation}
	\begin{split}
		& \hat{f}^t_{j_1,k_2} = \frac{1}{N} \sum\limits_{j_2=0}^{N-1} f(x_{j_1},x_{j_2};t) \exp{(-i k_2 x_{j_2})}  \\
		& \hat{f}^t_{k_1,k_2}= \frac{1}{N} \sum\limits_{j_1=0}^{N-1} f(x_{j_1},x_{k_2};t) \exp{(-i k_1 x_{j_1})}
	\end{split}
\end{equation}

可以看到，无论首先在哪一个方向上进行离散Fourier变换，最终得到的谱空间的Fourier系统都是相同的。
在进行并行时，从实空间到谱空间，或者从谱空间到实空间对矩阵进行的转置操作只需要进行一次。
这样可以减少数据在进程间的传递时间，从而提高计算速度。

\section{混淆误差}
在第二章第三节中，我们提到在谱空间中计算$\vec{u} \cdot \nabla \rho$和$\vec{u} \cdot \nabla \omega$的卷积消耗很大的计算量。
将其进行Fourier逆变换到实空间后，可以直接进行乘法运算得到其结果，然后再进行Fourier变换到谱空间去。
虽然通过这样的办法可以极大的降低计算量，但是这也很容易带来混淆误差。

事实上，在计算微分方程的非线性项时，使用上述方法都容易出现混淆误差。下面将简要介绍混淆误差的产生原因及处理方法。

\subsection{产生原因}
不失一般性，我们仅考虑一维情形，对$u$和$v$进行Fourier展开：
\begin{equation}
	u_j = \sum \limits_{k=-N/2}^{N/2-1} \hat{u}_k \exp(ikx_j), \quad v_j = \sum \limits_{k=-N/2}^{N/2-1} \hat{v}_k \exp(ikx_j), \quad j=0,1,\cdots,N-1.
\end{equation}
其中$x_j = 2 \pi j /N$，则$u$和$v$的Fourier系数为：
\begin{equation}
	\hat{u}_k = \frac{1}{N} \sum\limits_{j=0}^{N-1} u_j \exp(-ikx_j), \quad \hat{v}_k = \frac{1}{N} \sum\limits_{j=0}^{N-1} v_j \exp(-ikx_j) \quad k = -\frac{N}{2}, \cdots, \frac{N}{2}-1.
\end{equation}
并且定义：
\begin{equation}
	s_j = u_j v_j, \quad j = 0,1,\cdots,N-1.
\end{equation}
同样对$s$进行Fourier展开：
\begin{equation}
	s_j = \sum \limits_{k=-N/2}^{N/2-1} \tilde{s}_k \exp(ikx_j), \quad  j=0,1,\cdots,N-1.
\end{equation}
相应地Fourier系数为：
\begin{equation}
	\tilde{s}_k = \frac{1}{N} \sum\limits_{j=0}^{N-1} s_j \exp(-ikx_j),  \quad k = -\frac{N}{2}, \cdots, \frac{N}{2}-1.
\end{equation}

容易看出如下关系：
\begin{equation}
	\tilde{s}_k = \sum \limits_{m+n=k} \hat{u}_m \hat{v}_n + \sum \limits_{m+n=k\pm N} \hat{u}_m \hat{v}_n = 
	\hat{s} + \sum \limits_{m+n=k\pm N} \hat{u}_m \hat{v}_n
\end{equation}
其中等式右边的第二项即为混淆误差。

\subsection{处理方法——填充法}
在计算卷积项时如果有混淆误差的存在会造成解的不精确，处理混淆误差的方法主要有两种，填充法与平移法。
在本次程序中使用的是填充法。

填充法是指通过选取较多的网格点，在谱空间中将高频波的Fourier系数赋值为零，从而消除混淆误差。
下面将做更为具体的说明。

假设利用的网格点数为M，并且$M \ge N$，网格函数$\bar{u}_j$，$\bar{v}_j$进行Fourier展开：
\begin{equation}
	\bar{u}_j = \sum \limits_{k=-M/2}^{M/2-1} \breve{u}_k \exp(iky_i),  \quad
	\bar{v}_j = \sum \limits_{k=-M/2}^{M/2-1} \breve{v}_k \exp(iky_i),  \quad j = 0,1,\cdots,M-1.
\end{equation}
并且定义：
\begin{equation}
	y_j = 2\pi j/M  \quad
	\bar{s}_j = \bar{u}_j \bar{v}_j,  \quad j = 0,1,\cdots,M-1.
\end{equation}
这里，$\bar{u}$，$\bar{v}$，$\bar{s}$为$y_j = 2\pi j/M$的值，而$u_j$，$v_j$，$s_j$为$x_j = 2\pi j/M$的值。
其中，
\begin{equation}
	\breve{u}_k = 
	\left \{
	\begin{array}{ll}
		\hat{u}_k, \quad & |k| \le N/2 \\
		0, \quad & otherwise 
	\end{array}
	\right.
\end{equation}

那么，容易得到如下关系式：
\begin{equation}
	\breve{s}_k = \sum \limits_{m+n=k} \breve{u}_m \breve{v}_n + \sum \limits_{m+n=k\pm N} \breve{u}_m \breve{v}_n. 
\end{equation}
当$|M| > N/2$时，$\breve{u}_m = \breve{v}_n = 0$，所以等式右端的第二项为零。
这样就消除了混淆误差。

注意：为了保证$\breve{u}_m = \breve{v}_n = 0$，在选取网格点数时必须满足：
\begin{equation}
	-\frac{N}{2}-\frac{N}{2}\le \frac{N}{2}-1-M,\quad \text{即:}\quad M\ge\frac{3N}{2}-1
\end{equation}

\section{在时间方向利用高阶Runge-Kutta法离散布朗斯基方程}
在利用Fourier拟谱方法在空间上对布朗斯基方程进行离散后，得到了方程\ref{boussinesq_fft}。
这是一个关于时间t的常微分方程，需要在时间方向上继续进行离散。

由于在空间上的Fourier拟谱方法对Fourier基函数进行了截断，故损失了一定的精度，
这就需要我们在时间方向上使用更高精度的离散格式。
在这里，我们使用了现阶段最常用的求解常微分方程的方法——Runge-Kutta法。

\subsection{Runge-Kutta方法简介}
基于Taylor展开方法，从原则上可建立任意阶的单步差分方法，但是，由于它需要求出右端项的各阶偏导数形式，
所以在实际计算中这种直接使用 Taylor 展开来建立高阶差分格式的方法很少被采用。
为此，Runge 提出了间接使用 Taylor 展式来构造高精度数值方法的思想：
首先用函数f在m个点上的值的线性组合来代替y的导数，然后按Taylor级数展开，确定其中的组合系数。

这样既可避免计算y的高阶导数，又可保证差分方法具有较高的精度。
这使得Runge-Kutta法成为常用于求解常微分方程的一项技术，在计算中运用较为普遍。

\subsection{四级四阶古典显示Runge-Kutta公式}
按照构造法则可以构造出各种高阶Runge-Kutta方法，对于同级的Runge-Kutta 公式，隐式公式的精确度要比显示公式的高，
且又有较好的数值稳定性。但是在隐式Runge-Kutta公式中，由于需要解非线性方程组，这会给计算带来很大的开销和不便。
所以这里我们使用的是显式Runge-Kutta公式。

另一方面，由于Runge-Kutta方法的级数从4变为5时，精度的阶数不会理想地从4提高到5，
故4级以上的Runge-Kutta方法很少被采用。
这里采用的最常用的\textbf{四级四阶古典显式Runge-Kutta公式}：

假设有如下欧拉形式的常微分方程：
\begin{equation}
  \frac{d U}{d t} = L(t,U)
\end{equation}
其中U为待求函数，$L(t,U)$ 为右端项。
则相应的四阶的 Runge-Kutta 格式为: 
\begin{equation*}
	\begin{split}
   		&U_{n+1} = U_{n} +\frac{\triangle t}{6}(K_1+2K_2+2K_3+K_4), \\
  		&K_1 = L(t_n,U_n) \\
  		&K_2 = L(t_n+\frac{1}{2}\triangle t,U_n+\frac{1}{2}\triangle tK_1),\\
  		&K_3 = L(t_n+\frac{1}{2}\triangle t,U_n+\frac{1}{2}\triangle tK_2),\\
  		&K_4 = L(t_n+\triangle t,U_n+\triangle tK_3);
	\end{split}
\end{equation*}
令 $L(t,U) = \lambda U$，代入到上面的方程则：
\begin{equation*}
	\begin{split}
  		&K_1 =  \lambda U_n\\
  		&K_2 = \lambda (1+\frac{1}{2}\lambda \triangle t)U_n,\\
  		&K_3 = \lambda (1+\frac{1}{2}\lambda \triangle t+\frac{1}{4}\lambda^2 \triangle t^2)U_n\\
  		&K_4 = \lambda (1+\lambda \triangle t+\frac{1}{2}\lambda^2 \triangle t^2+\frac{1}{4}\lambda^3 \triangle t^3)U_n;
	\end{split}
\end{equation*}
从而有:
\begin{equation}
  	U_{n+1} = [1+\lambda \triangle t+\frac{1}{2}(\lambda \triangle t)^2+\frac{1}{6}(\lambda \triangle t)^3+\frac{1}{24}(\lambda \triangle t)^4]U_n
\end{equation}
故 4 阶 Runge-Kutta 方法的绝对稳定域为：
\begin{equation}
 	\mid 1+\lambda \triangle t+\frac{1}{2}(\lambda \triangle t)^2+\frac{1}{6}(\lambda \triangle t)^3+\frac{1}{24}(\lambda \triangle t \mid < 1
\end{equation}
事实上可以证明所有 m 阶 Runge-Kutta 方法的绝对稳定域都具有如上形式。

所以在实际计算中，为了满足稳定性，时间步长和空间步长必须满足相应的 CFL 条件，在本次程序中我是令$CFL = 0.5$。


\section{MPI并行程序思路}

\subsection{数据分块}
MPI并行算法必须考虑两个关键问题，即数据分块和数据通信。
合理的数据分块和数据通信将直接影响到MPI并行的性能。
本程序中使用的是一维条的数据分块，即对矩阵按行进行分块。
假设需要迭代的数据是一个$N\times N$的矩阵A(1:N,1:N)，进程的个数为$n\_rank$，
依次在每个进程上分配一个$M\times N,M = N/n\_rank$的数据块来存放数据，
且第i个进程上存放的数据为：$A(k+1:k+N,1:N),k=i\times M$。

\subsection{数据通信}
由于本程序使用的是FFTW一维串行函数来手工实现二维的 FFT，故对一个方向做完 Fourier 变换后，需要对另外一个方向进行 Fourier 变换，
此时需要对矩阵进行转置，并进行数据通信， 如下图所示：
\begin{figure}[!htp]
	\centering
	\includegraphics[width=\textwidth]{files/mpi.png}
	\caption{假设有两个进程，进程0和进程1。需要在x方向和y方向做Fourier变换。假设对x方向已经做完FFT，现需要在y方向做FFT，此时需要进行数据通信。其中左图为数据交换前的数据存储形式，右图为数据交换后的数据存储形式。}
\end{figure}

\begin{itemize}
\item 1. 绿色标定的数据块1在数据交换前后一直处在0号进程中，在数据交换后所有的数据的位置做了一次转置。
\item 2. 紫色标定的数据块4在数据交换前后一直处在1号进程中，在数据交换后所有的数据的位置也做了一次转置。
\item 3. 红色标定的数据块2在数据交换前处在0号进程中，在数据通信后传输到1号进程中，且对应数据的位置做了一次转置。 
\item 4. 蓝色标定的数据块3在数据交换前处在1号进程中，在数据通信后传输到0号进程中，且对应数据的位置也做了一次转置。 
\end{itemize}

在完成数据通讯后，只需要再在x方向做一次Fourier变换，即可完成两个方向的Fourier变换。

\subsection{混淆误差}
由前面的分析可以看出，在计算常微分方程的右端项时，需要用到逆离散 Fourier 变换先在实空间进行乘积计算， 
然后通过离散 Fourier 变换映射到谱空间。这中间会出现卷积的计算，容易出现混淆误差。
如果处理不当数值解会产生比较大的误差。

这里使用了填充法来处理误差，即在数组的中间插入一堆等于零的值，长度必须大于 1/2 的数组长度。
然而需要注意的是插零的位置——由于 FFTW 在处理 Fourier 变换时，利用了谱空间数组的共枙性质，将谱点在 $-N/2$
到 $-1$ 的值平移了一个周期（移到了 从 $N/2$ 到 $N$ 的位置），从而得到了一个下标为 $0,1,2,\cdots N-1$ 的数组。
所以零的位置必须放在数组的中间，而不是数组的两端。

注意：处理混淆误差可以把整个计算区域的剖分节点设为 3N/2，在进行 Fourier 变换时将数组中间 N/2 个数取为零，从而消去混淆误差。
但是我们这里不做这样的处理，因为在迭代过程中只有在计算右端項的时候需要处理混淆误差，需要在数组中间插入 N/2 个等于零的数，
这时数组的长度才会由 N 变为 3N/2。 但是在迭代的过程中其实不需要对数组这些插零位置的地方进行迭代。
为了减少迭代的计算量，每次计算完右端相后，再消去数组中间长度为 N/2 的等于零的数，把右端相的长度重新变回原长度 N 的大小。  

\subsection{Fourier变换的主要算法}
在利用常微分方程数值解法 Runge-Kutta 方法时，主要的运算量集中在利用当前步的 $\rho,\omega$ 的值来计算右端項的值，
这需要进行卷积的计算。 
如在计算 $\vec{u}\cdot\nabla\rho$ 和 $\vec{u}\cdot\nabla\omega$ 时，需要用到离散 Fourier 逆变换——先在实空间进行乘积计算，
然后通过离散 Fourier 变换映射到谱空间。这是整个程序的主要计算量所在。所以怎样实现一套高效的实空间的值和谱空间 Fourier 系数相互转换的算法， 对于程序的高效性至关重要。

下面给出了这两个过程的算法：

\begin{algorithm}
	\begin{spacing}{1.5}
	\caption{Fourier 变换 C2R}
	\begin{algorithmic}[1] %每行显示行号
		\Require \textbf{A}：存储谱空间的 \textbf{Fourier} 系数
    		\Ensure \textbf{B}：存储实空间的值
    		\Function {FFT\_C2R\_2D}{\textbf{A, B}}
    			\State 将 \textbf{A} 的值按行赋给一维数组 \textbf{in}.
    			\State 为了消除混淆误差, 在每个一维数组 \textbf{in} 的中间插入 N/2 长度的零.
    			\State 对每个一维数组 \textbf{in} 做谱空间到实空间的一维 \textbf{IFFT} 计算, 得到一维数组 \textbf{out}.
    			\State 用一个中间数组 \textbf{temp} 存储所有一维数组 \textbf{out} 的值.
    			\State 按 \textbf{5.2} 介绍的数据通信的方法, 对 \textbf{temp} 进行转置.
    			\State 将 \textbf{temp} 的值按行赋给一维数组 \textbf{in}.
    			\State 对每个一维数组 \textbf{in} 做谱空间到实空间的一维 \textbf{IFFT} 计算, 得到一维数组 \textbf{out}.
    			\State 用数组 \textbf{B} 存储所有一维数组 \textbf{out} 的值.
    		\EndFunction
  	\end{algorithmic}
	\end{spacing}
\end{algorithm}


\begin{algorithm}
	\begin{spacing}{1.5}
  	\caption{Fourier 变换 R2C}
  	\begin{algorithmic}[1] %每行显示行号
    		\Require \textbf{A}：存储实空间的值
    		\Ensure \textbf{B}：存储谱空间的 \textbf{Fourier} 系数
    		\Function {FFT\_R2C\_2D}{\textbf{A, B}}
    			\State 将 \textbf{A} 的值按行赋给一维数组 \textbf{in}.
    			\State 对每个一维数组 \textbf{in} 做实空间到谱空间的一维 \textbf{FFT} 计算, 得到一维数组 \textbf{out}.
    			\State 消去每个一维数组 \textbf{out} 因为混淆误差插入的网格点.
 	        	        \State 用一个中间数组 \textbf{temp} 存储所有一维数组 \textbf{out} 的值.
    			\State 按 \textbf{5.2} 介绍的数据通信的方法, 对 \textbf{temp} 进行转置.
    			\State 将 \textbf{temp} 的值按行赋给一维数组 \textbf{in}.
    			\State 对每个一维数组 \textbf{in} 做实空间到谱空间的一维 \textbf{FFT} 计算, 得到一维数组 \textbf{out}.
    			\State 消去每个一维数组 \textbf{out} 因为混淆误差插入的网格点.
   			\State 用数组 \textbf{B} 存储所有一维数组 \textbf{out} 的值.
    		\EndFunction
  	\end{algorithmic}
	\end{spacing}
\end{algorithm}
这样就完成了从谱空间的 Fourier 系数到实空间的值的运算。



\section{数值计算}

本文使用的并行环境是 MPICH-3.2。CPU 信息:Intel(R) Xeon(R) CPU E5-2650 @2.00 GHz × 16。
使用C语言来编写程序，使用FFTW包。

针对布朗斯基方程，使用周期性边界条件，在空间上使用FFT拟谱方法，在时间上使用四阶显式Runge-Kutta方法进行迭代求解。

\subsection{并行效率}
使用一维条进行数据分块，空间网格为$360\times 360$，时间尺度为$0.002$，运行到$t=2.5s$，总共迭代步数为$1250$步。
使用$1,2,4,8,16$个进程并行。

下表列出了使用不同进程数来运行程序的运行时间:
\begin{table}[htp]%有htp多个选择
\centering
\caption{运行时间}\label{tab:table} %用于ref
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
               \hline
	       进程个数 & 1 & 2 & 4  & 6 & 8 & 10 & 12 & 14 & 16  \\
	       %进程个数  &1 & 2  & 4  & 6 & 8 & 10 & 12 & 14 & 16 \\
               \hline
               %c & d \\
               运行时间(s) & 409.26 & 206.70 & 116.53 & 69.93 & 55.27 & 47.12 & 40.69 & 34.01 & 34.54 \\
               \hline
\end{tabular}
\end{table}

由于没有真解做比较, 我这里输出每个进程在2.5s 的数值解, 然后依次比较各组数据相差的误差, 结果达到了 1.0e-13, 基本达到了机器精度, 所以通过这样的方式来判断并行时进程间数据传递是否出错.

据此，可以计算出并行加速比和并行效率：

\begin{center}
加速比 = 串行执行时间 / 并行执行时间

并行效率 = 串行执行时间 / (并行执行时间 $\times$ 进程数)
\end{center}

\begin{figure}[h!]
	\center{
  	\includegraphics[width=0.48\textwidth]{files/accel_ratio.jpg}
	\includegraphics[width=0.48\textwidth]{files/accel_effects.jpg}}
  	\caption{并行程序的加速比与并行效率}
	\label{accel}
\end{figure}

由于我们的程序是在每个节点拥有16个核的Cluster上运行的，所以当$M \leq 16$时，
一方面，加速比都是随所用核数的增加而呈线性增加的。
另一方面，从加速效率图中不难看出，由于核与核之间的数据交换量较大，所以加速效率随着使用核数M的增加而降低。
除了在$M=4$处以及$M=16$处效率出现大幅下降之外，其余的加速效率基本稳定在$0.85$以上。

\subsection{数值结果}

下面画出了布朗斯基方程中密度 $\rho$ 在时间$t \in [0, 3.16]s$几个时间段的演化。空间计算区域为 $[0,2\pi]\times[0,2\pi]$，网格点个数为$512 \times 512$，使用周期边界条件，使用 Fourier 拟谱方法，时间方向迭代格式使用的是四级四阶Runge-Kutta法。

\begin{figure}[h!]
	\center{
  	\includegraphics[width=0.48\textwidth]{files/rho_0.png}
	\includegraphics[width=0.48\textwidth]{files/rho_700.png}
	\includegraphics[width=0.48\textwidth]{files/rho_1400.png}
	\includegraphics[width=0.48\textwidth]{files/rho_1900.png}}
  	\caption{数值结果$\rho$随时间的演化}
	\label{norm}
\end{figure}

\begin{figure}[h!]
	\center{
  	\includegraphics[width=0.9\textwidth]{files/rho_2000.png}}
  	\caption{数值结果$\rho$的发散结果}
	\label{abnorm}
\end{figure}

\begin{itemize}
\item 从图\ref{norm}中可以看出在 $t = 0s$ 时刻，初值为以 $y=\pi$ 为对称轴的对称图形，这与给定的初值有关。
在经过$t = 1.09s$到$t = 2.18s$后，图形像一个往右移动的水泡。
\item 随着时间的增加，图形沿着 x 方向移动，头部位置的形状变得更加尖锐，尾部的位置变得更加圆滑，且在尾部位置会留下一条又细又长的尾巴。
这个特征可以用来检验格式的耗散性，按理说低阶的耗散格式会破坏图形形状，不能明显的捕捉到这条尾巴。从物理上讲，有些类似KH不稳定性。
\item 但可以从图\ref{abnorm}中可以看出，当时间到$t=3.12$时，程序运算的结果开始发散，这时的解已经不能保证正确性了。
\end{itemize}

\newpage
\section{总结}
通过这次作业，我学习到了很多很多。

（1）因为采用的是显式计算格式，所以在程序编写方面省了不少时间。
然而由于必须满足CFL条件，在时间步长上存在较大的限制，如果网格划分的较精细，时间步长也必须足够小，因此计算需要花费很长时间。

（2）谱方法的要点是把解展开成光滑函数（一般是正交多项式）的有限级数展开式。
根据基函数的不同，有不同的展开方式。
这里由于是考虑周期性问题，故这里在空间上考虑以三角函数做为基函数的拟谱方法。
利用了一维的 FFTW 包来实现二维的 Fourier 变换。

（3）算法上的难点主要有三个：实空间的值和谱空间的 Fourier 系数的相互转化；各个进程间的数据通信；常微分数值格式的选取和计算。

（4）通过在不同平台上编译程序的经验，学习到了很多Linux系统下的知识。
比如配置PATH环境变量，安装软件并解决依赖树问题，查看CPU信息等等。

（5）在程序编译调试过程中遇到了很多问题，它们中有的给出了明确的出错地点，有的则并没有。
通过对这些问题的逐一排查并解决，学习到了编写并行程序时需要注意的很多细节，收获了很多经验。
程序的运行步骤参看包文件中的 README。

\begin{thebibliography}{99}
\bibitem{pre1} 李荣华，刘播，《微分方程数值解法(第四版)》，高等教育出版社，2009。
\bibitem{pre2} 张林波，迟学斌，莫则尧，李若，《并行计算导论》，清华大学出版社，2006。
\bibitem{pre3} 张平文，李铁军，《数值分析》，北京大学出版社，2007。
\end{thebibliography}









\end{document}


